Kinetic cross coupling between non-conserved and conserved fields in phase field 

models 



(N 

o 

(N 

S 1 

O 



o3 
o3 



O 
o 



m 
q 

o 

(N 



X 



Efim A. Brener 1 and G. Boussinot 1,2 
1 Peter Griinberg Institut, Forschungszentrum Jiilich, D-52425 Jiilich, Germany 
2 Computational Materials Design Department, Max-Planck Institut fur Eisenforschung, D-40074 Diisseldorf, Germany 

(Dated: July 3, 2012) 

We present a phase field model for isothermal transformations of two component alloys that in- 
cludes Onsager kinetic cross coupling between the non-conserved phase field cj> an d the conserved 
concentration field C. We also provide the reduction of the phase field model to the corresponding 
macroscopic description of the free boundary problem. The reduction is given in a general form. 
Additionally we use an explicit example of a phase field model and check that the reduced macro- 
scopic description, in the range of its applicability, is in excellent agreement with direct phase field 
simulations. The relevance of the newly introduced terms to solute trapping is also discussed. 



Introduction. Interface kinetics often plays a very im- 
portant role in phase transformations. It is responsible 
for the deviation from the local thermodynamic equilib- 
rium at the interface between different phases. In the 
case of alloys, it may affect the microstructure and the 
average concentration of the growing phase. It is respon- 
sible for solute trapping and some oscillatory instabilities 
of the transformation front leading eventually to the for- 
mation of banded structures. 

For example, in the case of isothermal transformations 
of binary alloys the bulk of each phase is described by the 
diffusion equation. The interface kinetics is more compli- 
cated because phase transformation effects occur in this 
region in addition to the diffusional exchange. In the 
macroscopic phenomenological approach to this problem 
(see, for example, [l[ and references therein) it is assumed 
that the physical interface width is much smaller than 
any relevant macroscopic length scale and one formulates 
linear Onsager relations at the interface which describe 
the interface kinetics. These relations connect two in- 
dependent fluxes J a and Jb (through the interface) of 
atoms A and B to two independent driving forces 5 ha 
and Shb which are the differences in the chemical po- 
tentials of A and B atoms at the interface (see below). 
The corresponding symmetric positive-definite Onsager 
matrix fully describes the interface kinetic properties in 
the framework of linear noncquilibrium thermodynam- 
ics. Of course, as in any phenomenological description, 
the three independent elements of this matrix depend 
on the specific physical mechanisms in the interface re- 
gion. This means that any specific thcrmodynamically 
consistent model of the interface (in general nonlinear) 
being linearized near the equilibrium must be reducible 
to this phenomenological description and the elements of 
the Onsager matrix should be calculated in terms of the 
model parameters. 

In recent years the phase field approach to phase trans- 
formations has attracted the attention of much research 
(see, for example, Q and references therein). It was orig- 
inally introduced as a mathematical tool to solve the free 
boundary problem without directly tracking the interface 



position. In the case of isothermal transformations of bi- 
nary alloys, this approach introduces, in addition to the 
conserved concentration field C, a nonconscrved phase 
field 4>. This field changes smoothly on the scale of the 
interface width from some value, say = 0, that corre- 
sponds to one phase to some other value, 0=1, which 
corresponds to the other phase. The phase field equa- 
tions of motion have a "diagonal" form in the classical 
variational formulation (see below), i.e. the time deriva- 
tive of <j> (C) depends only on the functional derivative 
of the free energy with respect to <j> (C). This "diagonal" 
formulation therefore contains only two independent co- 
efficients describing the interface kinetic properties while 
the general macroscopic phenomenology allows three in- 
dependent parameters. An intuitively clear way to re- 
solve this problem is to introduce kinetic cross coupling 
(non-diagonal terms) directly into the phase field equa- 
tions. However, to our best knowledge, this idea of a 
more general description of the interface kinetics in phase 
field models was never discussed in the literature (but see 
a very recent paper Q for a different approach). More- 
over, as stated in [2], according to Curie's principle [H, 
there can be no kinetic coupling between the scalar non- 
conserved phase field <f> and vectorial diffusional fluxes 
of the conserved quantities energy and/or concentration. 
We think that this is an erroneous statement (see also 
remarks in Q). The presence of the interface and the 
existence of the vector V0, that is orthogonal to the in- 
terface and operates only in the interface region, allows to 
formulate phase field equations that include kinetic cross 
coupling and are in agreement with linear nonequilibrium 
thermodynamics and Curie's principle. This issue is also 
relevant to the anti-trapping current introduced in some 
non- variational versions of the phase field model 0, HJ for 
different purposes. The anti-trapping current introduces 
a new kinetic coefficient and uses V0 as a vector nor- 
mal to the interface. To use this idea for the description 
of the cross effect of the interface kinetics in phase field 
models, one should carefully consider the necessary On- 
sager symmetry. This goal can be achieved only in the 
variational formulation of the phase field model. 
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The main purpose of this paper is to provide a ther- 
modynamically consistent description of the phase field 
model which contains cross kinetic coupling between the 
phase field <j> and the concentration field C. The second 
goal is to provide the reduction of this phase field model 
to the macroscopic phenomenological description as de- 
scribed above. This reduction is important to understand 
which macroscopic problem can be solved by the phase 
field model. The phase field model contains explicitly 
the finite interface width W as a parameter that is not 
included in the macroscopic description. We know only 
one example of such a reduction for classical ("diago- 
nal" ) phase field models which keeps W finite (the sharp 
interface limit sets W = 0). The thin interface limit 
was originally introduced by Karma and Rappel Q for 
the temperature field instead of the concentration field. 
This approach has then been promoted by many authors 
including the concentration field in the discussion (see, 
for example, a very detailed paper by Elder et. al [lCf). 
However, Brener and Temkin [6[ have recently pointed 
out that the macroscopic description derived by the thin 
interface limit has a clear deficit in some range of pa- 
rameters of the original phase field model. Namely, it 
can create strong unphysical instabilities due to a viola- 
tion of the positive-definiteness of the obtained Onsager 
matrix, while the original phase field model is fully con- 
sistent and stable. The approach promoted in the present 
paper for the reduction to the macroscopic description is 
free from this deficit. 

Macroscopic description of isothermal alloy transfor- 
mations. We discuss the phase transformation of a two 
component alloy at a given temperature T with an inter- 
face separating two phases. The dimensionlcss concentra- 
tion of B atoms is C\ in growing phase 1 and C% in mother 
phase 2. In the bulk of each phase these concentrations 
are described by diffusion equations with diffusion coef- 
ficients D\ and D 2 . In order to formulate the boundary 
conditions at the interface we use the phenomenological 
Onsager approach. Onsager relations connect the fluxes 
J A and Jb (through the interface) of atoms A and B to 
two driving forces Sfj, a and Sfis which are the usual differ- 
ences in the chemical potentials of A and B atoms at the 
interface (see, for example, [l| and references therein), 

6fiA/T = AJ A + BJ B , (1) 
6fi B /T = BJa+CJb ■ (2) 

The Onsager matrix should be positive-definite: A and 
C must be positive and B 2 < AC. According to the con- 
servation of B atoms at the interface we also have P, Hlj 

-Di(n- VCi) = VCx-Jb , (3) 
-D 2 (n ■ VC* 2 ) = VC 2 - Jb , (4) 
V = Ja + Jb , (5) 

where n is the unit vector normal to the interface and V 
is the normal velocity of the interface. In this description 



the matrix of Onsager coefficients describes a positive en- 
tropy production (per unit area), Ts — JaS^a + Jb5^b, 
in the interface region. For the following it is useful to use 
V = J a + Jb and Jb as independent fluxes and 5fj, A and 
5)i = Sf-iB — S/ia as corresponding driving forces. This 
choice preserves the invariance of the entropy production, 

Ts = Ja5^a + JbS^b = VS^la + JsSfi . (6) 

In this representation the linear relations between driving 
forces and fluxes read: 

Sfi A /T = AV + BJ B , (7) 
Sfi/T = BV+CJb , (8) 

with A = A,B = B + A,tmdC=C+A + 2B.We note 
that if f(C,T) is the free energy density of the phase, 
then often \i = /is — \ia = df/dC is called the diffu- 
sion chemical potential and [ia = /(C) — the grand 
potential. 

Phase field approach. We normalize the total free en- 
ergy F by T and write the dimensionlcss free energy G 
in the standard form for phase field models, 

G = F/T = J dv[H[W 2 {Vcj>) 2 /2+jWOW] +<?(<?,<«}. 

(9) 

Here fow = 2 (1~ 4 >2 ) is the normalized double-well po- 
tential which has equal minima at 4> = and <j> = 1; W 
is the characteristic scale of the interface width; g(C, <j>) 
is the dimensionlcss density of the chemical free energy; 
H represents the relative amplitude of the double-well 
potential normalized by T and H is usually a large pa- 
rameter. We assume that bulk phase 1 corresponds to 
d> = 1 and bulk phase 2 to </> = 0. Then, the dimension- 
less density of the bulk free energy is g(C,l) = f\{C)/T 
and g(C,0) = f 2 (C)/T. The detailed form of g{C,<f>) is 
model dependent. 

We write the system of phase field equations in the 
following variational form: 

- 5G/5(j) = T(j} + M$W(J ■ V0) , (10) 
-V{8G/5C) = M c W<j>V<f> + J/D((f>), (11) 
C+ (V • J) = . (12) 

In the bulk of each phase only the J terms survive lead- 
ing to the usual diffusional flux, Ji = — D{VCi and 
J2 = — U2VC2, with the bulk diffusion coefficients, D\ = 
[D{4> = l)/T)]dm/dC and D 2 = [D(cf> = 0)/T]dfx 2 /dC. 
In the interface region all terms are important leading to 
more complicated kinetics. 

The expression for the total entropy production reads 

S= J dV[-j>5G/8<j)-3 -V{5G/5C)} (13) 

= J dV[T(j>f + J 2 /D{(j)) + (Mef, + M C )W4>(3 ■ V0)] . 
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Onsagcr symmetry requires 



Ma = M c = M. 



(14) 



The conditions of positivc-definiteness of the entropy pro- 
duction read, r > 0, D(4>o) > and 



t/M 2 > max[W 2 D((l)o){V<f> ) 2 } , 



(15) 



where 0o is the phase field distribution at thermodynamic 
equilibrium. 

In classical phase field models M = 0. The terms with 
M represent kinetic cross coupling and introduce a new 
kinetic coefficient. A term analogous to our term with 
Mc has been introduced in 0, H| for a different purpose, 
using non-variational versions of phase field equations. 
To our best knowledge, the term with Ma has never been 
included before in phase field theory. Moreover, this term 
must be included in a thermodynamically consistent the- 
ory due to Onsager symmetry, Eq. (fl4j). as soon as the 
Mc term is included. 

Reduction of the phase field description to the macro- 
scopic description. We integrate the phase field equa- 
tions over the interface region in order to derive effective 
boundary conditions in the form of Eqs. (J7|) and (0. 
This will allow us to express the macroscopic elements 
of the Onsager matrix in terms of the phase field param- 
eters and to have an additional check of the symmetry 
condition. 

We assume that the interface is locally flat because 
we are mainly interested in kinetic effects rather than 
in the Gibbs-Thomson curvature correction and denote 
the direction normal to the interface by x. In the inter- 
face region we make a quasi-stationnary approximation, 
4> ~ —V(j>'(x) and C « — VC'(x), due to the strong gra- 
dients of and C in this region even at thermodynamic 
equilibrium. Integrating the continuity equation (TT2"]) in 
the interface region and choosing the integration constant 
equal to —Jb we find, 



J(x) w - J B + VC{x). 



(16) 



Eq. (|16|) then reproduces macroscopic continuity equa- 
tions, (J3j) and (|4]), if the observation point x is chosen in 
phase 1 or 2 near the interface. We integrate Eq. (fTT|l . 



M C W / dx\^{x)] 2 - / dx C (x)/D((t> Q ) 
Jw Jw 



5n/T = V 

+ J B I dx/D{<f) ) 



(17) 



w 



and also, multiplying Eq. (|10p by (j>'(x) and integrating 
over the same region, we find 

Sfi A /T = V\r [ dx[cf>' (x)} 2 + f dx C 2 (x)/D(M 
1 Jw Jw 



(Mf + M c )W dx[^{x)] 2 C {x) 



w 



+Jb 



MaW dx[fi (x)} 2 - dxC (x)/D(fo) 
Jw Jw 



(18) 



Here J w denotes the integral over the interface region 
whose width is of order W, but such that <j) ranges from 
4> ~ 1 to (f> « 0. We have replaced cf>(x) and C(x) by 
their equilibrium distributions, <fio(x) and Cq(x), due to 
linearization. We have also used the following steps in or- 
der to integrate the left- hand-side of Eq. (ITU|) . First, the 
contribution proportional to H vanishes (l2|. Secondly, 
we write 



dx 



4>'{x)dg/d<j)= [ dg- [ dx C \x)dg/dC , 
Jw Jw 



<W JW JW 

integrate the last term by parts 



dx cj>'{x)dg/d(t> = S^a/T + dx C(x){dg/dC)' , 



w 



and use again Eq. (fTTj) for (dg/dC)'. 

As expected Eq. (fTBj) has the form of Eq. and Eq. 
([T7]) has the form of Eq. ((5J . The Onsager symmetry of 
Eqs. (0 and (jHJ requires Ma = M c = M that confirms 
Eq. (fT4"]) . We can also easily check that the interfacial 
part of the entropy production, Eq. (|13[) . reduces to the 
form of Eq. © with 5fi given by Eq. (|17p and 8 fi a given 
by Eq. (fT8|) . It is of course positive if the condition (fT5"|) 
is fulfilled. The phase field model presented here con- 
tains three independent inverse velocity scales describing 
the interface kinetics, t/W, J w dx/D((f>) and M, while 
classical phase field models include only two. 

Explicit example and numerical checks. Our aim now 
is to compare quantitatively simulation results within 
a specific phase-field model to the solution of the cor- 
responding macroscopic description using the reduction 
presented above. We focus on the one-dimensional 
steady-state growth of phase 1 at the expense of phase 
2, a case where the growth velocity V is kinctically con- 
trolled. Due to the global conservation law, the concen- 
tration C\ in phase 1 is constant, C\ — Coo, where C^o 
is the concentration far ahead of the interface in phase 2. 
Within the macroscopic description, in the limit of small 
velocity, V and Ci read 0] 



V 

c 2 



(/l'(Cr 9 )/T)(C 1 e9 -C 00 )AC 



BiC\ q 

+ (Coo 



) + cci q c e 2 q 

Cl q ) + {B + CGl q )V 



(19) 
(20) 



where /{'(C) is the second derivative of fi(C) with re- 
spect to C and AC* = C 2 q - C\ q with C eq (C e 2 q ) the two- 
phase equilibrium concentration of phase 1(2). 

We use a simple phase-field model for which the chemi- 
cal free energy densities fi{C) and /2(C) of phases 1 and 
2 parabolically depend on the concentration, 



g(C,<b) = ±(c-C e 2 q ~ P ( ( j))Ac) 



(21) 



with p((f>) = 3 (1O - 150 + 6</> 2 ) (see for example jl3|). 
For an equilibrium interface centered at x — 0, we have: 
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<t> (x) = 1/2 - taah[x/(y/2W)]/2 with O = 1 in phase 
1 and O = in phase 2; C (z) = {C eq + C e /)/2 + 
u(x)AC/2 with u(x) = -u{-x) = 1 - 2p[<p (x)}. 

For simplicity, and in order to make further analyti- 
cal progress, we assume a constant diffusion coefficient 
D((j)) = D. This assumption is physically more relevant 
to solid-solid transformations than to solidification prob- 
lems where D\ <C D 2 . We perform the integrations in 
Eqs. (fTTl ITS)) in a symmetric range [—5, 6] around x = 
yielding: 



A 

B 
C 



w -f3WAC 2 /(4D) 

+ [(Cf) 2 + {C e 2 q ) 2 ]5/D - Ma(CZ q 
Ma-(C? + C c 2 q )5/D , 
25/ D, 



d q 



)(22) 
(23) 
(24) 



where 5 ~ W but such that cf>o(—S) « 1 and </>o(<$) ~ 0. 
Then 

da; [0 o (x)] 2 ~ WM cfc [</> <»] 2 ~ 0.23570 , 

-5 J — oo 

P = f g ^ [1 - « 2 (z)] « |" ^ [1 - U 2 (z)] « 1.40748 , 

due to the fast convergence of the integrals. Using the 
latter kinetic coefficients in Eqs. (|19[) and (|20[) . we find 
FandC 2 (/i'(Cf)/T=l): 



V = 



(CT* - C oc )AC 



ar/W - f3WAC 2 /(4D) ' (25) 
C 2 = C 2 " + Coo - + (M a - SAC/D) V. (26) 

While the velocity is essentially independent of the in- 
tegration range i5 as discussed in the concentration 
C 2 depends on 5. We note that in our reduced descrip- 
tion the interface concentrations and chemical potentials 
are actually defined at the spatial points x = ±5 and 
vary slightly with 5 due to weak gradients if the system 
slightly deviates from equilibrium. For a more detailed 
discussion of this issue and its relation to the extrapola- 
tion procedure in the thin interface limit see 0. 

In Fig. [TJ we present a comparison of the dimension- 
less velocity Vt/W as a function of given by the 
analytical formula, Eq. (|2"5"j) . and obtained from phase- 
field simulations. The two equilibrium concentrations are 
Cl q = 0.3 and C 2 q = 0.7, the diffusion coefficient (con- 
stant throughout the whole system) is Dt/W 2 = 0.5 and 
H = 50 (we checked that the results are essentially in- 
dependent of H for such large values). We find a good 
quantitative agreement in the linear regime, i.e. for small 
velocities. The simulations reproduce the independence 
of M for the velocity in the linear regime (see Eq. (|25l) ) . 
Nonlincaritics of the phase field model naturally lead to 
deviations at higher velocities. We mention that here the 
denominator in Eq. (|25[) is positive. For smaller values 
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FIG. 1: Dimensionless velocity Vt/W vs concentration of the 
system Coo for different values of M (crosses: MW/t — 2; 
circles: MW/t = 0) compared with the analytical prediction 
(line) of Eq. (J5SJ. The case MW/t = -2 is indistinguishable 
from MW/t = 2. 
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FIG. 2: Partition coefficient k vs dimensionless velocity for 
different values of M (crosses: MW/t = 2; circles: MW/t = 
0; triangles: MW/t = — 2) compared with the corresponding 
analytical prediction (lines) of Eq. ()26|) with 8 — 2y/2W. 



of D, the denominator may be negative and steady-state 
solutions exist even for Coo > C[ q (see, for example, 
0, EH EH and references therein) . 

In Fig. [2] we present the partition coefficient, k = 
Ci/C 2 (we recall that in steady-state C\ = Coo), as a 
function of the dimensionless velocity for different values 
of M, with C 2 measured at x = 6 = 2\[2W . The classi- 
cal phase field model (M = 0) shows already the solute 
trapping effect (increase of the partition coefficient with 
velocity) while positive values of M show anti-trapping 
tendency and negative values of M promote further so- 
lute trapping. This was the reason for authors of 0, H| 
to include a term with positive Mq in Eq. (|ll[) calling 
it anti-trapping current. However, we understand now 
that a thcrmodynamically consistent description requires 
simultaneously to include the term with — Mq in 
phase field Eq. ([TO]). 

We have also checked numerically the stability con- 
dition, Eq. (|15|) (for our explicit example it reads 
8t/(DM 2 ) > 1), by investigating the relaxation to the 
equilibrium configuration. If the condition is violated by 
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1.5% the system "blows up" instead of relaxing to the 
equilibrium. 

Summary. We have formulated a phase field model 
given by Eqs. (JTUJ), (fTTj) . and (fT2|) . It includes Onsager 
kinetic cross coupling between the non-conserved phase 
field 4> and the conserved concentration field C. We have 
performed the reduction of this model to the correspond- 
ing macroscopic description given by Eqs. (fTT]) and (fl"8|) . 
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